#ifndef DIRKN_MATRIX_HPP
#define DIRKN_MATRIX_HPP

#include <vector>
#include <algorithm>
#include <numeric>
#include <stdexcept>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <assert.h>

namespace tws {

    template <typename T>
    class Diagonal
    {
        protected:
            const std::size_t rows_;
            const std::size_t cols_;
            std::vector<T> data_;
        public:
            typedef T value_type;
            Diagonal(std::size_t rows, std::size_t cols, const T& init)
            : rows_(rows)
            , cols_(cols)
            , data_(std::min(rows, cols))
            {
                std::fill(data_.begin(), data_.end(), init);
            }
            const T/*&*/ operator()(std::size_t i, std::size_t j) const {
                return (i == j) ? data_[i] : 0; // can't return reference to const 0
            }
            const T& operator()(std::size_t i) const {
                return data_[i];
            }
            T& operator()(std::size_t i) {
                return data_[i];
            }
    };
    
	template <typename T>
	class ColumnMajor
	{
		protected:
			const std::size_t rows_;
			const std::size_t cols_;
			std::vector<T> data_;
		public:
			typedef T value_type;
			ColumnMajor(std::size_t rows, std::size_t cols, const T& init)
			: rows_(rows)
			, cols_(cols)
			, data_(rows*cols)
			{
				std::fill(data_.begin(), data_.end(), init);
			}
			const T& operator()(std::size_t i, std::size_t j) const{
				return data_[j*rows_+i];
			}
			
			T& operator()(std::size_t i, std::size_t j) {
				return data_[j*rows_+i];
			}
	};
	
	template <typename T>
	class RowMajor
	{
		protected:
			const std::size_t rows_;
			const std::size_t cols_;
			std::vector<T> data_;
		public:
			typedef T value_type;
			RowMajor(std::size_t rows, std::size_t cols, const T& init)
			: rows_(rows)
			, cols_(cols)
			, data_(rows*cols)
			{
				std::fill(data_.begin(), data_.end(), init);
			}
			const T& operator()(std::size_t i, std::size_t j) const{
				return data_[i*cols_+j];
			}
			
			T& operator()(std::size_t i, std::size_t j) {
				return data_[i*cols_+j];
			}
	};
	template <typename T>
	class Symmetric
	{
		protected:
			const std::size_t rows_;
			const std::size_t cols_;
			std::vector<T> data_;
		public:
			typedef T value_type;
			Symmetric(std::size_t rows, std::size_t cols, const T& init)
			: rows_(rows)
			, cols_(cols)
			, data_((rows+rows*cols)/2)
			{
				assert(rows==cols);
				std::fill(data_.begin(), data_.end(), init);
			}
			const T& operator()(std::size_t i, std::size_t j) const{
				int k=abs(j-i);
				return data_[(-1)*k*(k-1-2*cols_)/2+std::min(j,i)];
			}
			
			T& operator()(std::size_t i, std::size_t j) {
				int k=abs(j-i);
				return data_[(-1)*k*(k-1-2*cols_)/2+std::min(j,i)];

			}
	};

    template <class Shape>
    // Op deze plaats gebeurt de overerving!!, matrix erft alles van een template superklasse!
    class Matrix : public Shape
    {
        public:
            typedef typename Shape::value_type value_type;
            Matrix(std::size_t rows, std::size_t cols, const value_type& fillwith = 0) 
            : Shape(rows, cols, fillwith)
            {
            }
            // Deze methodes zijn nodig om rows_ en cols_ te kunnen opvragen
            std::size_t rows() const { return this->rows_; }
            std::size_t cols() const { return this->cols_; }

            static int print_width;
            friend std::ostream& operator<<(std::ostream& out, const Matrix& A)
            {
            		std::cout << std::setprecision(std::numeric_limits<long double>::digits10+1)
					<< std::scientific;
                for(std::size_t i=0; i < A.rows(); ++i) {
                    for(std::size_t j=0; j < A.cols(); ++j)
                        out << std::setw(print_width) << A(i,j) << " ";
                    out << std::endl;
                }
                return out;
            }
    
            
            template <class MatrixConcept>
            Matrix& operator=(const MatrixConcept& A)
            {
            	assert((this->rows_ == A.rows()) && (this->cols_ == A.cols())); 
            	for(std::size_t i=0; i < A.rows(); ++i) {
            		for(std::size_t j=0; j < A.cols(); ++j) {
						(*this)(i,j) = A(i,j);
            		}
        		}
        		return *this;
            }
    };
    
	
    template <class Shape>
    int Matrix<Shape>::print_width = 30;


    
    template <class Matrix>
    typename Matrix::value_type max_norm(const Matrix& A)
    {
        using std::abs; using std::max;
        auto s = A(1, 1);
		for(std::size_t i=0; i < A.rows(); ++i) {
            for(std::size_t j=0; j < A.cols(); ++j) {
				s = max(abs(A(i,j)),s);
            }
        }
        return s;
    }
            		
    

    
    template <typename T>
    T max_norm(const Matrix<ColumnMajor<T>>& A)
    {
        using std::abs; using std::max;
        auto s = A(1, 1);
		for(std::size_t j=0; j < A.rows(); ++j) {
            for(std::size_t i=0; i < A.cols(); ++i) {
				s = max(abs(A(i,j)),s);
            }
        }
        return s;
    }

}

#endif
